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Abstract. We introduce and experimentally demonstrate a technique for performing 

quantum state tomography on multiple-qubit states despite incomplete knowledge 
about the unitary operations used to change the measurement basis. Given unitary 
operations with unknown rotation angles, our method can be used to reconstruct 
the density matrix of the state up to local az rotations as well as recover the 
magnitude of the unknown rotation angle. We demonstrate high-fidelity self-calibrating 
tomography on polarization-encoded one- and two-photon states. The unknown unitary 
operations are realized in two ways: using a birefringent polymer sheet — an inexpensive 
smartphone screen protector or alternatively a liquid crystal wave plate with a 
tuneable retardance. We explore how our technique may be adapted for quantum 
state tomography of systems such as biological molecules where the magnitude and 
orientation of the transition dipole moment is not known with high accuracy. 
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1. Introduction 

Quantum state characterization is essential to the development of quantum technologies, 
such as quantum computing [HE], quantum information |3] and quantum cryptography 
[3]. The measurement of multiple copies of the quantum state and subsequent 
reconstruction of the state's density matrix is known as quantum state tomography 
(QST) [51 i. 

QST relies on performing measurements in a number of different bases to gain 
access to all the information required to reconstruct the density matrix. Different 
measurement bases are typically obtained by applying unitary operations to the state 
before measurement. For systems that couple to the electromagnetic field via their 
transition dipole, precise knowledge of the magnitude and orientation of the dipole 
moment is vital to performing such measurement-basis changes. For example, QST of 
trapped ion qubits requires precise control of the orientation of the ions via magnetic 
fields [7]. 

These techniques can not be immediately applied to systems that are either not well 
characterized, such as biological molecules, e.g. those found in photosynthetic systems 
[El |9], where the magnitude and orientation of the transition dipole moment is not 
known with high accuracy; or systems that experience variability during the fabrication 
process, such as colloidal quantum dots [TO] . 

With such systems in mind we introduce a method for self- calibrating tomography 
(SCT), which successfully reconstructs the state of the system despite incomplete 
knowledge about the unitary operations used to change the measurement basis. Other 
efforts to characterise photosynthetic systems are also underway, e.g. quantum process 
tomography using ultrafast spectroscopy was proposed by Yuen- Zhou et al. [TH [T2] and 
Hamiltonian tomography was proposed by Maruyama et al. [13j . 

Our method for self-calibrating tomography may also improve the robustness of 
standard tomography, where calibrated unitary operations undergo small errors and 
fluctuations over the course of the experiment. For other modifications to standard 
QST due to inaccessible information or preferable measurements choices, we refer the 
reader to references [a[ISl[ini[IZl[IHllini[2ni[nil221[2Sll2a[2Sl[2SlE^ 

In Section |2| we show that given local unitary operations with unknown rotation 
angles, but known and adjustable rotation axes, it is possible to reconstruct the density 
matrix of a state up to local rotations, as well as recover the magnitude of the 
unknown rotation angles. 

We demonstrate SCT in a linear-optical system using polarized photons as qubits 
in Section |3] An inexpensive smartphone screen protector, i.e. an uncharacterized 
birefringent polymer sheet, is used to change the measurement basis. We go on to 
investigate the technique's robustness to measurement noise and retardance magnitude, 
and demonstrate SCT of a two-qubit state using liquid crystal wave plates with tuneable 
retardances. 

In section [4], we investigate the application of this technique to quantum state 
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tomography of photosynthetic systems. 
2. Theory 

In this section, we give a brief review of quantum state tomography and introduce our 
technique for self-cahbrating tomography. 

2.1. Quantum State Tomography 

The state of a qubit, given by the density matrix p, can be decomposed into a sum of 
orthogonal operators 0",: 

1 ^ 

P = 2 5Z ^^^^ ' (1) 

i=0 

where (Tq is the identity operator and ai_3 are the Pauh matrices. The coefficients Xj 
are given by the expectation values of the orthogonal operators, Aj = Tr[/3(Tj]. 

Given an unknown state, p can be reconstructed by performing a number of 
projective measurements fij on subsequent copies of the state. The measurement 
statistics are given by 



n 



Ar,Tr[p/i,] , (2) 

where Mj is a constant that depends on the duration of data collection, detector 
efficiency, loss etc. Inserting ([T]) into ^ gives 

^^■ = TE^^Tr[/i,^,]. (3) 

i=0 

The equations given by ([s]) can be solved for Aj, which can then be inserted into ([T]) 
to reconstruct the unknown density matrix p. The operators fij must be chosen such 
that the equations in (|3]) form a linearly independent set. If Mj are known, three 
different measurements suffice, given the normalization condition Aq = 1. Typically, 
Mj is unknown, but assumed to be equivalent for all measurements, in which case, four 
measurements are required to reconstruct the state. In practice, due to fluctuations in 
measurement statistics, it is advantageous to employ a maximum likelihood estimation 
method [5] to reconstruct p rather than directly solving (|2]), which can lead to unphysical 
density matrices. 

We can deflne the measurement operators (ij in terms of a fixed known measurement 
operator (xq and different unitary operations Uj as follows. 

Pj = UjpoUj . (4) 

In terms of the measurement statistics, this is equivalent to applying the unitary operator 
Uj to the state prior to making measurements given by Pq: 



n 



AfjTiipUjpoU,] = Af,TT[UjPUjM , (5) 



where we used the cyclic properties of the trace. 
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2.2. Self- calibrating Tomography 

Now consider the case where the unitary operator, used to change the measurement 
basis, is characterized by an unknown parameter a. Recall that a unitary operation is 
simply a rotation on the Bloch sphere. In this paper, we restrict ourselves to rotations 
about axes that lie on the equatorial plane of the Bloch sphere, i.e. the x—y plane, 
given by 

If. (a) = (ua) = Q-i'^»i^os{^)-^+sm{v)ay)/2 

= cos ( ^ ) (To - 2 sin ( ^ ) {cos{ip)cr^ + sm{if)ay) , (7) 



V 2 / ^ V 2 

where the index j refers to a particular (f and z/, where ip is the known angle between 
the axis of rotation and the and u is real number that, along with a, determines 

the rotation angle. As we see later, the x—y plane restriction arrises naturally in the 
systems we consider. The measurement operators are defined as 

jj,j{a) = Uj{a)fioUj{a) . (8) 

The relationship between the measurement statistics and the parameters which 
characterize the density matrix will be given by 

nj = ^J2^T[fi,{a)a,]\. (9) 

i=0 

By making an additional measurement, it is possible to solve the series of equations given 
by ([9]) for Aj and \a\. As before, measurements are chosen such that generate a 

linearly independent set of equations. If this condition is met, the equations yield two 
sets of solutions, corresponding to ±a, i.e. either the actual state or the phase-flipped 
version of the state. 

It is straightforward, but non-trivial, to extend this formalism to multiple-qubit 
systems. However this is not required as in practice one may use a maximum likelihood 
estimation method adapted for SCT to reconstruct the density matrix, which we will 
introduce in the next section. 



2.3. Maximum Likelihood 

The solutions for Aj in ^ and ^ are typically quite sensitive to fluctuations in 
measurement statistics, which can lead to unphysical density matrices. It is therefore 
favourable to use a maximum likelihood algorithm to reconstruct the density matrix. 

Assuming the noise on the measurements has a Gaussian probability distribution, 
the probability of obtaining a set of N measurement results {uj} is [5] 

N _(nj-n^)2 

P = ylJJe '^'^ , (10) 

j=0 

where A is a normalization constant and aj is the standard deviation of the 
jth measurement (given approximately by y/nj). The expectation value of each 
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measurmeent is given by 

=^iTr[p/ij(a)] , (11) 

where p is the state we are trying to solve for. To ensure a normahzed, positive and 
hermitian density matrix, we restrict p to the form 



where 



P 



T 



Tr[ftf] 
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(13) 



for a single qubit and 
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(14) 



for a two qubit system. Maximising the likelihood that the physical density matrix 
p gave rise to the data {uj} — i.e. maximising (10) — can be reduced to finding the 
minimum of the function 

^ ' - ^2 ^(n,-Ar,Tr[pA,(a)])^ 

j=0 



£({t,},a) 



j=0 



[rii 



2a] 



2n^ 



(15) 



For large uncertainty in the measurement statistics, the likelihood function can 
have multiple local minima and the maximum likelihood algorithm may converge to a 
state and retardance that do not correspond to the actual state. This can be overcome 
given a priori knowledge of the purity of the state, or by taking more measurements. 



3. Linear-optical implementation 

For a linear-optical test of this method, we performed SCT on polarization-encoded one- 
and two-qubit states. 



3.1. Single qubits 

To prepare single-qubit states, we first pumped a 1mm long type-I down-conversion 
beta-barium borate (BBO) crystal with a 405nm diode laser, as shown in Fig. [T] 
a), which leads to spontaneous generation of pairs of photons. The detection of a 
horizontally polarized photon in one spatial mode of the down-converted photon-pair 
heralds the presence of a horizontally polarized photon in the other mode, i.e. the single 
qubit. Using the method described in Bocquillon et al. [31], and based on an estimated 
average rate of 2-photon events in the signal conditional on a detected photon in the 
heralding arm of 0.96s~^, we calculated the g^'^\0) of our heralded single-photon source 
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Figure 1. Experimental scheme: a) A single-qubit state is prepared, unknown unitary 
operations are implemented using BPSs or an LCWP and the state is measured using 
a trusted, fixed measurement; b) An entangled two-qubit state is prepared, unknown 
local unitary operations are implemented using LCWPs and the two-qubit state is 
measured using a trusted, fixed measurement. 



to be ^''•^''(O) = 0.067, indicating a high quahty single-photon source. Quarter- and 
half-wave plates can then prepare the qubit in arbitrary polarization states. 

For the trusted, fixed projector, we chose fio = \R) {R\, implemented with a quarter- 
wave plate and a polarizing beamsplitter, followed by coincidence- counting, as shown 
in Fig. [l] a). To change the measurement basis, a unitary operation was applied to 
the state prior to measurement, implemented using either one or two identical pieces 
of a birefringent polymer sheet (BPS). The BPS is an optical element that alters the 
polarization state of light in the same way as a wave-plate, with the crucial difference 
that the retardance of the BPS is not known. The effect of the BPS on the state of the 
photon can be described by R(f,{a) given in ([6|, where a is the retardance of the BPS, 
is the alignment of the BPS's optical axis with respect to the polarization of the light, 
and o"j are in the \R) /\L) basis. Using different alignments of the BPS, we constructed 
the following operators /ij(a;) = tjj{ay fioOj^a) using different implementations of tjj{a): 



f/o(a) = ao, (16) 

Ui{a) = Ro{a) , (17) 

f/2(«) = ^|(«), (18) 

f/3(«) = /?.(2«), (19) 

Ui{a) = R3.{2a) . (20) 



Note that rotations of 2a were required to satisfy the linear-independence condition 
on These rotations were realized with a double-pass through the BPS — the sheet 
was cut in half and the pieces placed in succession. The beam passed through sections 
of the pieces that were within close proximity before partition, to ensure consistency in 
the birefringence. One piece was then removed to perform rotations of a, before removal 
of the remaining piece to measure /to- This measurement sequence eliminated the need 
for realignment between measurements. 

In this particular linear-optical implementation of the scheme, both output ports 



Self- calibrating Quantum State Tomography 



7 




Figure 2. Single-shot fidelity between SCT and ST for 14 states distributed on the 
Bloch sphere, using: a) a birefringent polymer sheet (BPS) as the uncalibrated wave 
plate; and b) using a liquid crystal wave plate (LCWP) set to different values of the 
retardance a. c) Histogram showing the predicted retardance and the fidelity between 
SCT and ST, for 100 runs of SCT with varying numbers of total photons collected 
(and hence levels of counting error). All used \4>) — \H) and the BPS. 



of the PBS are easily accessible. We therefore measured the complimentary projectors, 
fi'j{a) = Uj{ay fiQUj^a) where fiQ = \L) (L|, to determine Mj. 

Using this method, SCT was performed on 14 states: six that lie at the ends 
of the three axes and eight other states defined by all combinations of 6 = ±| and 
= {±|,±^} where \ip) = cos (|) \R) + e*''^sin (|) \L). Each state was reconstructed 
using the maximum-likelihood estimation method described in Section 2^ in conjunction 
with the YALMIP solver [32]. 

For comparison, standard tomography (ST) was also performed on each of the 
states. ST required fewer projectors fij = UjfioUj, given by the unitary operators: 

(21) 
(22) 



Uo 


= CTo 










Ui 


= Ro{ 












= Rk 


(D 




2 





(23) 

and their complimentary projectors, fi'j = UjfiQUj where fiQ = \L){L\. We used the 
fidelity F = Ti[>y v^PstPsctVpst]^ to characterize the success of SCT, where psT 
and psCT are the density operators reconstructed using ST and SCT respectively. We 
achieved an average fidelity F = 0.99 ± 0.01 and determined the retardance of the BPS 
(at 405nm) to be a = 0.58 ± 0.03 (close to 7r/6). Individual fidelities for each state are 
shown in Fig. |2| a) . 
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3.2. Dependence on retardance 

To determine the effectiveness of SCT as a function of tlie retardance, we repeated the 
experiment using Bolder Vision Optik hquid crystal wave plates (LCWP), as shown in 
Fig. [l]a). In an LCWP, the birefringence is a function of the applied voltage, allowing 
for the retardance to be set arbitrarily. A single LCWP was used to implement rotations 
of a and 2a using different voltage settings. Fidelities for different setting of a are shown 
in Fig. [2]b). The average fidelity decreases with the size of the retardance. The fidelity 
for states lip) = \R) and I?/') = |L) always remains close to unity. This is a consequence 
of our choice of measurement projectors. Using the LCWP, we find a slightly lower 
average fidelity with a greater spread over the different states, even for values of a 
larger than that of the BPS. We attribute this to errors associated with the non-linear 
dependence of the retardance on the voltage passed through the LCWP. This influences 
how accurately 2a can be set with respect to a. 

3.3. Dependence on noise 

To investigate the dependence of the fidelity on the amount of noise, we performed SCT 
on the state {ip) = \H) using a retardance a = 7r/32. The noise scales inversely with 
the square root of the number of counts and by collecting from 150 to 3000 photons we 
varied the noise between roughly 1% and 10%. For each noise setting, SCT was repeated 
100 times. The resulting fidelities and predicted retardances are presented in Fig. |2]c). 
Notice that the majority of predicted fidelities and retardances correspond to the actual 
values. However, recall that the likelihood function can have multiple local minima, 
making it possible for the maximum likelihood algorithm to converge to a state and 
retardance that do not correspond to the actual state. This shows up in the high-noise, 
or low photon-number, region and can be seen as a small bump in the fidelity histogram 
(and to a lesser extent in the retardance histogram) in Fig. [2]c). As the noise decreases, 
only high-fidelity results corresponding to the actual state and retardance remain. Fig. 
[3] shows the distribution of states on the Bloch sphere as predicted by SCT for high and 
low noise amounts. The darker red cluster of states in Fig. |3] a) also correspond to a 
local minimum in the likelihood function. 

3.4. Entangled qubits 

For multi-qubit tomography, entangled two-qubit states, lip) = a\HH) + h\VV) , were 
prepared using two 1mm long type-I down-conversion BBO crystals with optical axes 
aligned in perpendicular planes [33] pumped by a 405nm diode laser, as shown in Fig. 
[l]b). The parameters a and h which characterise the state were tuned by changing 
the polarization of the pump beam with a HWP. LCWPs in each arm, implemented 
the unknown unitary operations, followed by detection of (iq = |-R) (-R| in both modes. 
For this two-qubit experiment, LCWPs were chosen over the BPSs due to realignment 
issues that accompany the insertion (as opposed to removal) of optical elements in the 
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Figure 3. Distribution of states on the Bloch sphere as predicted by SCT for: a) 
high noise; and b) low noise. Notice the darker red cluster of states in Fig. a), which 
corresponds to a local minimum in the likelihood function. There are 100 states shown 
in each sphere. 



beam-path. The concurrence [3l] for states reconstructed using SCT and ST, as well as 
the fidelity between the two, are shown for a variety of entangled states in Table [T] 



SCT 



a 



ST 



0.905 ±0.005 
0.566 ±0.005 
0.337 ±0.005 
0.004 ±0.004 



0.927 ±0.003 
0.562 ± 0.005 
0.328 ± 0.004 
0.010 ±0.003 



0.968 ± 0.0028 * 
0.978 ±0.0012 
0.977 ±0.0032 
0.989 ±0.0015 ** 



Table 1. The concurrence [34] as predicted by both tomographic techniques 
and fidelity between the two reconstructions for states with different degrees of 
entanglement. Density matrices for * and ** are shown in Figure |4] a) and b) 
respectively. Error bars were determined by using Poissonian noise in a Monte Carlo 
simulation 

Fig. |4] shows density matrices for the most- and least-entangled states. The 
retardances were determined to be ai = 0.810 ±0.009 and a2 = 0.760 ±0.006, compared 
with the known value a = tt/A ^ 0.785. 



4. Potential applications to photosynthetic systems 

In this section, we give some detail about how this technique may be applied to 
performing quantum state tomography on biological molecules, such as chromophores 
in photosynthetic systems, where the transition dipole moment and its orientation are 
not known with high accuracy. 



4.1. Unitary operation 

Consider a molecule with two states of interest — the ground state \g) and the excited 
state |e) — that undergoes an interaction with an EM field, e.g. a laser pulse, given by 



Self- calibrating Quantum State Tomography 10 




Figure 4. Comparison between density matrices reconstructed using SCT and ST: a) 
F = 0.968 ±0.0028; h) F = 0.989 ± 0.0015. For concurrences and other fidelities, refer 
to Table m 



the interaction Hamiltonian 

H = d- E{t) = 4 ■ E{t)\e) {g\ + dge ■ E{ty\g) (e| , (24) 

where E{t) is the electric field, d is the transition dipole moment operator and 
deg = {e\d\g) = d*^. An ultra-fast pulse with a carrier frequency Uc can be described 
by E{t) = S{t)e^'^''^e + h.c, where S{t) is the temporal profile of the pulse and e is 
the polarization vector. Under the rotating wave approximation, we can write the 
Hamiltonian in the interaction picture as 

H = (4 ■ e)£(t)e*(--'^'=)*|e) {g\ + (4 ■ e}*£itre-'^'^^-^^^'\g) (e| , (25) 

where Ue is the \g) — ?■ |e) transition energy. The unitary governing the interaction is 
given by 

Uj{a) = e"^/'^*^^ = g-»(Re[n]CT^+Im[0]aj;)/2 ^ ^-iua{cos{ip)a^+sm(v)ay)/2 ^ j^26) 

where Q = l{deg ■ e) J dtS{t)e'^'^--'^^^^ and 
2 

" = 11 

a = (4 ■ e) , (28) 
= arg ( / dtS{t)e'^'^^-'^^A . (29) 



dt£{t)e 



i{u}c — U)e)t 



(27) 



Comparing with (pi), we can see that (26) constitutes a rotation of angle z/a, where 



(y9 — determined by the relative phase of £(t) — defines the rotation axis. Notice from 



(27) that V is linearly proportional to the square root of the intensity of the pulse at the 
transition frequency. We can therefore execute rotations of 2a by changing the intensity 
of the pulse, without any knowledge of deg ■ e. 

4-2. Projective measurement 

We can perform a measurement corresponding to the operator 

= |e) (e| , (30) 
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by preparing multiple copies of the state and performing statistics on the number of 
photons emitted. A molecule in the excited state will — on a long enough time-scale — 
emit a photon, while one in the ground state will not; a molecule in a superposition 
of ground and excited states a\g) + b\e) will emit a photon with probability The 
measurement statistics will be given by 

Uj = f/jTT[Uj{a)pUj{ay fio] , (31) 

where Uj{(y) is given in (26) and Afj depends on the solid angle of the detector, the 
position of the detector with respect to the dipole moment of the molecule and losses. 

Measurement statistics can be collected either by repeated preparation and 
measurement of a single molecule, recently demonstrated by Hildner et al. [321 ES], 
which is analogous to the linear-optical implementation discussed above; or by preparing 
an ensemble of molecules in the same state and measuring the fluorescence of the 
sample, where the time-integrated intensity of the emitted light is directly proportional 
to the number of photons. An open question is whether it will be possible for the self- 
calibrating technique to compensate for an ensemble of randomly oriented molecules. 
This is currently being investigated. 

In principle, the extension to multiple entangled molecules should be straightfor- 
ward: even if individual molecules can not be spatially isolated, independent rotations 
on each qubit could be executed with appropriately shaped pulses, as long as there 
are different transition energies for each molecule. However, in the case of coupled 
molecules, there are technical difficulties associated with addressing the transition ener- 
gies of the individual molecules — the system rather absorbs light at the eigen-energies of 
the coupled Hamiltonian. In this system of two coupled two-level systems can be 

approximated as a three- level system (where the fourth level is typically disregarded). 
Self-calibrating quantum state tomography of coupled two-level system will require an 
extension of our formalism to a three-dimensional Hilbert space. This is also currently 
being investigated. 



5. Conclusion 



In conclusion, we have introduced and demonstrated a technique for performing 
quantum state tomography of multi-qubit systems which does not rely on complete 
knowledge of the unitary operation used to change the measurement basis. Our 
technique characterizes the state of interest with high fidelity as well as recovering the 
unknown unitary operation. We find that in the context of QST of polarization-encoded 
qubits, it is possible to do away with a well-calibrated half-wave plate and replace it 
with an inexpensive uncharacterized piece of birefringent material. 

A future extension of this work will attempt to incorporate unknown rotation axes 
and non-unitary transformations, such as those that introduce noise to the system, 
thereby reducing the purity of the state by an unknown amount. 

We anticipate our method to be applicable to quantum state tomography of 
systems that do not have well-characterized transition dipole moments, such as biological 
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molecules like those found in photosynthetic systems, as well as systems that experience 
variability during the fabrication process, such as colloidal quantum dots. 

AMB, AD and DFVJ thank DARPA (QuBE) for support. DHM, LAR and AMS 
thank NSERC, CIFAR and Quantum Works for support. We also thank Xingxing Xing 
for the serendipitous discovery of the birefringence of his smartphone screen protector. 
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